rm(list = ls())
library(sf)
library(tidyverse)
library(terra)
library(nhdplusTools)
library(mapview)
library(dataRetrieval)
library(lubridate)
library(prism)
library(ggspatial)
library(nngeo)# Added from OG code
library(stars)# Added from OG code
# this gives you an error, but it can be ignored:
try(plyr::ldply(list.files(path="src/",
pattern="*.R",
full.names=TRUE),
source))
FALSE Error in list_to_dataframe(res, attr(.data, "split_labels"), .id, id_as_factor) :
FALSE Results must be all atomic, or all data frames
# Rmarkdown options
knitr::opts_chunk$set(echo = T, warning = F, comment = F, message = F)
# mapview options
mapviewOptions(basemaps.color.shuffle=FALSE,basemaps='OpenTopoMap')
For this code to run properly, your site data must be configured as follows:
site.longitude
and latitude. Knowledge of coordinate projection required.
OR: Each site has their known COMID, with column name
comid.data/
folder.I have included an example data set called
placeholder.csv, along with all of the additional data sets
necessary for the code to run, stored in the data
folder.
Currently, this workflow requires downloading several data sets
locally for much speedier run times. This includes: PRISM climate &
aridity rasters, NHD flow direction data, and CONUS-wide NHD catchments.
All data sets are found in the shared data folder.
This analysis is only appropriate for locations along adequately-sized streams. Some streams are too small to be captured by NHDPlusV2; it is also common for coordinates to fall in the wrong catchment (especially for big rivers). For that reason, review each site and make sure that the NHD feature attributes were appropriately captured.
Identify each sample’s NHD COMID. This COMID will allow site linkages
with all datasets in this workflow. If COMID is already listed in the
CSV, make site_type = "comid".
sf_use_s2(FALSE)
site_type = "xy" # OR site_type = "comid"
sites <- read_csv("data/WHONDRS_WROL.csv")%>%
distinct(site, .keep_all = TRUE) %>%
dplyr::select('site','latitude','longitude')
Additional steps for sites with coordinates, and no COMID:
if(site_type == "xy"){
sites <- sites %>%
dplyr::select(site, latitude, longitude) %>%
sf::st_as_sf(coords = c("longitude","latitude"), crs = 4269) # 4269 = NAD83 CRS
if(sf::st_crs(sites) != sf::st_crs(4269)){
sites <- sites %>% st_transform(., crs = 4269)
}
mapview(sites)
}
Pull all meta data associated with each site’s COMID.
if(site_type == "xy"){
sites <- getNHDxy(df = sites)
}
FALSE [1] "45 locations linked to the NHD."
if(site_type == "comid"){
sites <- getNHDcomid(df = dplyr::select(sites, site, comid))
}
Make NHD-based watershed shapefiles for all CONUS sites. To make this
step MUCH faster, it is best to have a locally downloaded version on the
National NHD catchment shapefile stored on your local system. I have
already included this shapefile in the data folder.
site_watersheds <- getWatersheds(df = sites, make_pretty = TRUE) %>%
inner_join(., dplyr::select(sf::st_drop_geometry(sites), site, comid), by = "comid")
Ensure that each site is accurately represented by the NHD (some
locations may be located on streams that are too small to be captured by
the NHD, others may have coordinates that are slightly off, potentially
placing them in the wrong catchment). Here, we create simple maps of
each watershed, to determine if the delineated watershed/NHD attributes
are appropriate. This can take a while depending on how many sites you
are doing this for. Maps are automatically stored in the
data/maps/ folder. It is highly recommended to review each
map, particularly for known locations along BIG rivers and TINY
streams.
map2(sites$site, sites$comid, getMaps)
FALSE [[1]]
FALSE NULL
FALSE
FALSE [[2]]
FALSE NULL
FALSE
FALSE [[3]]
FALSE NULL
FALSE
FALSE [[4]]
FALSE NULL
FALSE
FALSE [[5]]
FALSE NULL
FALSE
FALSE [[6]]
FALSE NULL
FALSE
FALSE [[7]]
FALSE NULL
FALSE
FALSE [[8]]
FALSE NULL
FALSE
FALSE [[9]]
FALSE NULL
FALSE
FALSE [[10]]
FALSE NULL
FALSE
FALSE [[11]]
FALSE NULL
FALSE
FALSE [[12]]
FALSE NULL
FALSE
FALSE [[13]]
FALSE NULL
FALSE
FALSE [[14]]
FALSE NULL
FALSE
FALSE [[15]]
FALSE NULL
FALSE
FALSE [[16]]
FALSE NULL
FALSE
FALSE [[17]]
FALSE NULL
FALSE
FALSE [[18]]
FALSE NULL
FALSE
FALSE [[19]]
FALSE NULL
FALSE
FALSE [[20]]
FALSE NULL
FALSE
FALSE [[21]]
FALSE NULL
FALSE
FALSE [[22]]
FALSE NULL
FALSE
FALSE [[23]]
FALSE NULL
FALSE
FALSE [[24]]
FALSE NULL
FALSE
FALSE [[25]]
FALSE NULL
FALSE
FALSE [[26]]
FALSE NULL
FALSE
FALSE [[27]]
FALSE NULL
FALSE
FALSE [[28]]
FALSE NULL
FALSE
FALSE [[29]]
FALSE NULL
FALSE
FALSE [[30]]
FALSE NULL
FALSE
FALSE [[31]]
FALSE NULL
FALSE
FALSE [[32]]
FALSE NULL
FALSE
FALSE [[33]]
FALSE NULL
FALSE
FALSE [[34]]
FALSE NULL
FALSE
FALSE [[35]]
FALSE NULL
FALSE
FALSE [[36]]
FALSE NULL
FALSE
FALSE [[37]]
FALSE NULL
FALSE
FALSE [[38]]
FALSE NULL
FALSE
FALSE [[39]]
FALSE NULL
FALSE
FALSE [[40]]
FALSE NULL
FALSE
FALSE [[41]]
FALSE NULL
FALSE
FALSE [[42]]
FALSE NULL
FALSE
FALSE [[43]]
FALSE NULL
FALSE
FALSE [[44]]
FALSE NULL
FALSE
FALSE [[45]]
FALSE NULL
Interactive map showing all sites and their delineated watersheds:
mapview(site_watersheds, col.regions = "#56B4E9", alpha.regions = 0.2, lwd = 3, layer.name = "Watershed") +
mapview(sites, cex = 5, col.regions = "black", layer.name = "Points") +
mapview(st_read('data/site_flowlines.gpkg', quiet = T), lwd = 3, color = "red", layer.name = "Flowline")
If after reviewing the delineated watersheds you have found that some site coordinates place that location in the wrong catchment, we will need to explore that site’s nearby catchments to link it to the correct COMID:
# HJA-WS1, HJA-WS2, LOO-BLU
sus_points <- sites %>%
# list weird sites here:
filter(site %in% c("HJA-WS1","HJA-WS2","LOO-BLU"))
sus_nhd <- sus_points %>%
sf::st_buffer(., dist = 0.01) %>%
sus_mapper(.)
# Add a point for the USGS Gauge station
station = as.data.frame(46.53429376)
colnames(station) = "Latitude"
station$Longitude = -120.4672855
sbux_sf <- st_as_sf(station, coords = c("Longitude", "Latitude"), crs = 4326)
mapview(sus_nhd[[2]], col.regions = "#56B4E9", alpha.regions = 0.4, legend = FALSE) +
# mapview(sbux_sf, cex = 5, col.regions = "red", legend = FALSE) +
mapview(sus_points, cex = 5, col.regions = "black", legend = FALSE) +
mapview(sus_nhd[[1]], lwd = 1.5, color = "red", legend = FALSE)
# EAS-ROC, EAS-RUS, EAS-BRA, EAS-PUM
sus_points <- sites %>%
# list weird sites here:
filter(site %in% c("EAS-ROC", "EAS-RUS", "EAS-BRA", "EAS-PUM"))
sus_nhd <- sus_points %>%
sf::st_buffer(., dist = 0.01) %>%
sus_mapper(.)
# Add a point for the USGS Gauge station
station = as.data.frame(46.53429376)
colnames(station) = "Latitude"
station$Longitude = -120.4672855
sbux_sf <- st_as_sf(station, coords = c("Longitude", "Latitude"), crs = 4326)
mapview(sus_nhd[[2]], col.regions = "#56B4E9", alpha.regions = 0.4, legend = FALSE) +
# mapview(sbux_sf, cex = 5, col.regions = "red", legend = FALSE) +
mapview(sus_points, cex = 5, col.regions = "black", legend = FALSE) +
mapview(sus_nhd[[1]], lwd = 1.5, color = "red", legend = FALSE)
# CT-W9 and CT-POPE
sus_points <- sites %>%
# list weird sites here:
filter(site %in% c("CT-W9","CT-POPE"))
sus_nhd <- sus_points %>%
sf::st_buffer(., dist = 0.01) %>%
sus_mapper(.)
# Add a point for the USGS Gauge station
station = as.data.frame(46.53429376)
colnames(station) = "Latitude"
station$Longitude = -120.4672855
sbux_sf <- st_as_sf(station, coords = c("Longitude", "Latitude"), crs = 4326)
mapview(sus_nhd[[2]], col.regions = "#56B4E9", alpha.regions = 0.4, legend = FALSE) +
# mapview(sbux_sf, cex = 5, col.regions = "red", legend = FALSE) +
mapview(sus_points, cex = 5, col.regions = "black", legend = FALSE) +
mapview(sus_nhd[[1]], lwd = 1.5, color = "red", legend = FALSE)
# CT-BUNN too small no COMID?
sus_points <- sites %>%
# list weird sites here:
filter(site %in% c("CT-BUNN"))
sus_nhd <- sus_points %>%
sf::st_buffer(., dist = 0.01) %>%
sus_mapper(.)
# Add a point for the USGS Gauge station
station = as.data.frame(46.53429376)
colnames(station) = "Latitude"
station$Longitude = -120.4672855
sbux_sf <- st_as_sf(station, coords = c("Longitude", "Latitude"), crs = 4326)
mapview(sus_nhd[[2]], col.regions = "#56B4E9", alpha.regions = 0.4, legend = FALSE) +
# mapview(sbux_sf, cex = 5, col.regions = "red", legend = FALSE) +
mapview(sus_points, cex = 5, col.regions = "black", legend = FALSE) +
mapview(sus_nhd[[1]], lwd = 1.5, color = "red", legend = FALSE)
# CT-PHEL and CT-NEPA have same COMID. not accurate
sus_points <- sites %>%
# list weird sites here:
filter(site %in% c("CT-PHEL","CT-NEPA"))
sus_nhd <- sus_points %>%
sf::st_buffer(., dist = 0.01) %>%
sus_mapper(.)
# Add a point for the USGS Gauge station
station = as.data.frame(46.53429376)
colnames(station) = "Latitude"
station$Longitude = -120.4672855
sbux_sf <- st_as_sf(station, coords = c("Longitude", "Latitude"), crs = 4326)
mapview(sus_nhd[[2]], col.regions = "#56B4E9", alpha.regions = 0.4, legend = FALSE) +
# mapview(sbux_sf, cex = 5, col.regions = "red", legend = FALSE) +
mapview(sus_points, cex = 5, col.regions = "black", legend = FALSE) +
mapview(sus_nhd[[1]], lwd = 1.5, color = "red", legend = FALSE)
# # Based on the map, it appears that this site should be linked to comid a different comid
# updated_sites <- tibble(site = c("site8"),
# comid = c(17416032))
#
# sites <- updated_sites %>%
# getNHDcomid(.) %>%
# # Replace this new data in the sites df:
# bind_rows(filter(sites, !site %in% sus_points$site))
#
# site_watersheds <- getWatersheds(df = updated_sites, make_pretty = TRUE) %>%
# inner_join(., select(sf::st_drop_geometry(updated_sites), site, comid), by = "comid") %>%
# bind_rows(filter(site_watersheds, !site %in% sus_points$site))
If after reviewing the delineated watersheds you have found that some locations are in fact not appropriate for this analysis, remove them now:
# small_sites <- c("site9")
#
# # site 9 appears to be a headwater stream too small for analysis with NHDPlus V2:
# sites <- sites %>%
# filter(!site %in% small_sites)
#
# site_watersheds <- site_watersheds %>%
# filter(!site %in% small_sites)